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Abstract — Recursive feedback is defined and discussed 
as a framework for development of specific algorithms 
and procedures that propagate the time-domain solution 
for a dynamical system simulation consisting of multiple 
numerically coupled, self-contained, stand-alone subsys- 
tem simulations. A satellite motion example containing 
three subsystems (orbit dynamics, attitude dynamics, 
and aerodynamics) has been defined and constructed 
using this approach. Conventional solution methods are 
used in the subsystem simulations. Centralized and dis- 
tributed versions of coupling structure have been ad- 
dressed. Numerical results are evaluated by direct com- 
parison with a standard total-system, simultaneous- 
solution approach. 

L INTRODUCTION 

Digital simulations of dynamical systems are often built 
by constructing algorithms that solve a set of differential- 
algebraic equations that mathematically model the system. 
The equations are solved numerically, as a single coupled 
set, using one of several standard or modified numerical 
integration methods. In some cases, software is written in a 
language such as FORTRAN or C to define the equations, 
and an existing or slightly modified ordinary differential 
equation solver, perhaps a Runge-Kutta implementation, is 
employed to perform the integration. In other cases, an in- 
house development such as Marshall System for Aerospace 
Simulation (MARSYAS) or a commercial off-the-shelf 

product such as MATLAB^/SIMULINK® is used to act as 
a higher level facilitator of what is ultimately a mathemati- 
cally similar approach. 

System-level digital solution of coupled “stand-alone” 
dynamical simulations is a fundamentally different approach 
to system simulation, and fundamentally different numerical 
procedures are required. Recursive feedback is a conceptual . 
method from which a family of appropriate algorithms, 
processes, and specific numerical procedures can be de- 
rived. 

Coupled subsystem simulation architecture provides 
near complete independence of stand-alone subsimulations, 
and naturally facilitates high fidelity and broad scope 


through collaboration across interfaces that can be imple- 
mented in the same physical and engineering terms that de- 
fine them in the actual system. Such simplicity promotes 
clarity of communication and ease of understanding, both of 
which have many positive benefits. Individual sub- 
simulations can potentially be implemented on separate, 
remote, and dissimilar computational platforms, and this 
portends numerous advantages and possibilities as computer 
network capabilities improve. 

H. RECURSIVE FEEDBACK ALGORITHM 
A. Concept 

Referring to Fig. 1, the subsystems are represented by 
independent, self-contained, time-domain, dynamical simu- 
lations that can be commanded to run from a given set of 
initial conditions over a predetermined segment of time 
when provided with input signals as functions of time over 
that interval. In general, the time segment, also referred to 
as a convergence interval, is short compared with the total 
duration of the simulation, but may be significantly larger 
than the integration step size associated with a typical sub- 
system. For each new stage of recursion, revised subsystem 
input signals are computed by summing system-level input 
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Figure 1: Schematic of Recursive Feedback Process 







signals with current stage feedback (coupling) signals that 
have been generated as if the loop were open. New subsys- 
tem responses and feedback signals are computed using the 
revised subsystem input signals in conjunction with original 
initial conditions. The process is started by ignoring the 
feedback signal and continues recursively until convergence 
is achieved. After convergence, initial conditions are re- 
placed with final conditions and the process is restarted for 
the next segment of time. System response over the total 
desired duration of the simulation is the concatenated set of 
combined subsystem responses over many time segments. 

B . Insights for Specific Design 


where the parenthesized superscript is the stage index. Next, 
the stage zero feedback function is computed by multiplying 
the stage zero output response function by k. After “mov- 
ing” the feedback function numerical data across the recur- 
sive data transfer point (future stage to cuiTent stage in Fig. 
1) and summing it with the system-level input function 
(zero), it becomes the stage one subsystem input function. 
Stage one response can now be computed by numerical in- 
tegration of the new subsystem input function with applica- 
tion of the original initial condition y (0) . In the absence of 
numerical error, symbolic representation of the result is 

J>> (1) « = )’(0)[l + fe] . (4) 


Because Fig. 1 describes a system of considerable gen- 
erality whose complexity for a specific case can range from 
simple to great, further clarification of the process seems 
appropriate. The simplest “system” could arguably be a 
single integrator in the forward path, a constant multiplier 
(gain) in the backward (feedback) path, no (zero) system- 
level input, and a nonzero initial condition on the output. In 
that case, there is one linear subsystem, and that subsystem 
is a single integrator. The output function is the integral of 
the subsystem-level input function with respect to time and 
can be generated independently and directly by simply 
computing area under the curve and adding the initial condi- 
tion. The coupling matrix becomes a scalar (value of the 
loop gain) and actually represents direct feedback rather 
than coupling. The process, however, does not explicitly 
make such a distinction, nor does it need to. The summing 
junction becomes moot since there is no system-level input. 
The point at which feedback data are transferred from one 
stage to the next remains as depicted in Fig. 1 but does not 
exist in the analog system to be digitally simulated. The 
scalar differential equation that models the analog system is 


Repeating the process, symbolic representation of the stage 
two result is 


;(2) (0 = ;y(0)[l +£? + ~r~\ . 


By induction, the nth stage response is given by 
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Clearly, if sufficiently continued, the symbolic recursive 
process yields a solution that approaches the exact solution 
(2). In theory, the time interval (convergence interval) for 
this case can be as long as desired; in practice, it must be 
restricted because of numerical integration and round-off 
errors. 


As another example, consider the “system” to be an un- 
damped oscillator governed by the scalar differential equa- 
tion 


y(0 + © 2 y(0 = 0, (7) 


y(0 = ty (0 . (i) 

where t is time and k is the loop gain, and the exact closed- 
form solution for the output response is given by the expo- 
nential function 

y(*)=y(0)[«"] , (2) 

where y(0) represents the initial value of the output y(t). 

Applying recursive feedback to this system, the output is 
first computed as a function of time over the convergence 
interval as if there were no (zero) feedback. This function is 
referred to as the stage zero output response. Because the 
system input function is zero for this example, the response 
is simply a constant function of time symbolically repre- 
sented by 


where CO is the frequency of oscillation. A simulation of 
this system by recursive feedback can be set up in at least 
two very different ways. One design possibility is to define 
a single subsystem, best described as a double integrator, or 
two single integrators in series. The output function of the 
subsystem is determined by integrating its input function 
twice. The coupling matrix once again degenerates to a sca- 
2 

lar constant, — CD , and the system input function is zero. 
Assuming a nonzero initial condition y(0) on the output 
(displacement), application of recursive feedback to this 
system results in a truncated series definition of y(f) , and 
at the nth recursion stage, a symbolic representation is 
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Then for a sufficiently large number of recursion stages, the 
approximate solution approaches the exact solution 

y(t) = y(0)[COSQ)t] . (9) 


Once again the number of terms in the series corresponds 
directly to the number of recursion stages, but the order of 
the series with respect to time is twice the number of recur- 
sion stages. 

In contrast, a second design possibility is to define two 
subsystems that are both single integrators. The governing 
differential equation is now given by the two-dimensional 
matrix-vector equation 


m=KY(t), do) 


where the coupling matrix is given by 


K = 



and the output function is given by 
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where (t) represents displacement and y 2 (0 represents 
velocity. Application of recursive feedback to this system 
results in a truncated series definition of Y ( t ) . At the nth 
recursion stage, the approximate system response function 
is symbolically represented by 


V <n) (t) = 


Again, each recursion stage, in effect, adds one term to the 
truncated series approximation of the exact solution. 

For the two-subsystem, twin-integrator approach, it 
takes two recursion stages to reach the same order of ap- 
proximation in time that was reached in one stage with the 
single-subsystem, double-integrator approach. However, the 
twin integrators represent independent operations that can 
be numerically performed in parallel, while the double inte- 
grator represents either a single operation, or two operations 
that must be performed in series. In the double integrator 
case, velocity, y(f) , does not appear at the system level, 

only one function crosses the recursive data transfer point, 
and only one function is available to be checked for conver- 
gence at the system level. In the twin integrator case, veloc- 
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ity, ^ 2 ( 0 * appears at the system level, y^and 

y 2 (0both cross the recursive transfer point, and both 

functions are available to be checked for convergence at the 
system level. Clearly, in a more complex situation, the rela- 
tive merits of these approaches become simulation design 
and performance issues. 

For the broader range of systems contained within the 
framework of Fig. 1, the nature of the coupling matrix is a 
function of many possible simulation design choices. In its 
most basic form, a coupler (element of the coupling matrix) 
would simply route signals from one subsystem to another 
without modifying them. In its most complex form, it could 
become a multidimensional nonlinear operator as well. As 
implied by the double integrator example when viewed as 
two single-integrator subsystems in series, it is also possible 
to distribute coupling functions throughout the system. 
Thus, the concept of a coupling matrix is not essential to 
the process; it may, however, be useful in organization and 
control of a large simulation at the system level. In a par- 
ticular situation, one might choose to view several subsys- 
tem simulations as a single combined set when the nature of 
the coupling among them is such that an output response 
can be determined directly from an input signal without 
necessity for recursion among members of the set. Typical 
coupler functions are likely to be coordinate transforma- 
tions, gain multiplications, interpolations, and other data 
modification or routing tasks that must be performed to 
make proper interface connections among a set of prede- 
fined or existing subsimulations. 

C. Convergence 

To be of practical value in simulation of dynamical sys- 
tems, the ability of the recursion process to achieve conver- 
gence needs to be understood in a general sense so that rea- 
sonable assurance of convergence can be provided for spe- 
cific circumstances. The convergence issue has been ad- 
dressed for some simple nonlinear examples as well as a 
time- varying example in [1]. For the nonlinear examples, 
each recursion stage adds higher order terms; but for a 
given number of correct terms in the series expansion, a 
greater number of recursion stages is required. For these 
systems, the convergence interval cannot be arbitrarily 
long, even without integration and round-off errors, be- 
cause the series becomes divergent. The degree of restric- 
tion required for convergence depends on the initial condi- 
tion. 

The analysis of [2] addresses linearly coupled subsys- 
tems where the subsystems are also linear but are multidi- 
mensional. While subsystem output signals are linearly 
related to the subsystem states, the number of output signals 


Aerodynamic Data: 

Drag Coefficient 225 
Density 1.95 x 10 “ n kg/m 3 


may be less than the number of states. The recursive feed- 
back process is analytically applied based only on inputs 
and outputs to the subsystems, as depicted in Fig. 1. Be- 
cause no algebraic feed-through of subsystem input signals 
is allowed, the possibility of system-level algebraic loops is 
excluded. However, the analysis of [3], though much less 
straightforward, addresses a similar system with inclusion 
of system-level algebraic loops. 

From this work and other experience, it appears that 
convergence can be achieved for a broad range of ordinary 
systems simply by controlling the length of the conver- 
gence interval. It is apparent, however, that this is not al- 
ways the case. For example, systems with lower gain alge- 
braic loops may converge, while those with higher gain 
may not. The essence of this problem is illustrated by sub- 
stituting a unit gain multiplication process for the integra- 
tion process in the example system described by (1). This 
creates a purely algebraic loop. After adding a nonzero in- 
put function, analytical application of the recursive feed- 
back algorithm generates a geometric series rather than an 
exponential series, and convergence becomes conditional. It 
should be noted that a nonconvergent system of this type 
might also be physically untenable or not meaningful. The 
"‘guarantees” and insights with respect to convergence that 
come with analyses like that in [3] are weakened at best and 
totally lost at worst without the linear model assumption. 
However, one does have a guarantee of solution existence 
and uniqueness for a broad range of sufficiently well be- 
haved nonlinear systems [4], and a converged recursive 
feedback process is indicative of a solution with exception 
of questions concerning discretization error, round-off er- 
ror, inappropriate tolerances, etc. By shortening the length 
of the convergence interval, the number of recursion stages 
required to achieve convergence is normally decreased; or, 
in a nonconvergent situation, the likelihood of convergence 
is increased. The possibility exists for automatic in-process 
control of the convergence interval length as well as the 
number of recursion stages, not unlike control of stepsize 
and choice of order in a conventional algorithm for numeri- 
cal solution of differential equations. 

m. EXAMPLE SYSTEM 

The example system consists of a satellite in low-Earth 
orbit that is experiencing aerodynamic effects in addition to 
gravitational effects. The spacecraft is idealized as a rigid 
body of rectangular “box” shape, illustrated in Fig. 2. The 
center of mass is offset from the geometric center by a small 
amount in all three axes, and moments of inertia are such 
that there are two large but somewhat unequal principal 
values, while the third principal value is much smaller. The 
principal axes are normal to the box surfaces so that all 
cross products of inertia are zero. The models are based on 



width. 4 m 
height: 3 m 


Cwr Qf Mass Offset- 

X* 1.0m 
Y = 0. 1 m 
Z = 0 1 m 


Notes: 

1. 0 depicts center of mass location. 

Z rail, pitch, and yaw are Euler angles about the 
X, Y, and Z -axes, respectively. The sequence 
is yaw -pitch 'roll 

3 /, w, and A represent length, width, and height. 

4. All axes penetrate facets at their geometric centers. 


Ma ts Properties : 

Mass 10,000 kg 
Moments of Inertia- 
'll ,250.0 0 0 

0 193,125.0 0 

0 0 1,731,8125. 


kg “m 2 


Figure 2: Spacecraft Properties and Coordinate System 


constant density atmosphere, “panel” aerodynamics, and 
spherical Earth gravity. There is no attitude control, so the 
body is allowed to tumble under the influence of gravity and 
aerodynamics. Gravitational and aerodynamic torques affect 
the rotational dynamics of the rigid body, while gravita- 
tional and aerodynamic forces affect the translational dy- 
namics. Both translational and rotational dynamics affect 
aerodynamic forces and torques, so a natural (analog) feed- 
back loop is apparent. This example falls in the general 
category of nonlinear, continuous systems. 

IV. INDIVIDUAL SUBSYSTEM SIMULATIONS 


A. Orbit Dynamics 


The orbit dynamics simulation input signal is aerody- 
namic force, and output signals are orbit radius and velocity 
vectors. All orbit dynamics simulation variables are refer- 
enced to an Earth-centered inertial coordinate system as 
defined in Fig. 3. The underlying mathematical model, of 
the type found in [5], is represented by 
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where G is the universal gravitational constant, 
and m sat are masses of the Earth and satellite, and R is 

the orbit radius vector. R extends from the center of the 
Earth to the center of mass of the satellite. F aero is the 

aerodynamic force, and is derived from output of the aero- 
dynamics subsystem simulation. After rearrangement to 
first-order form (six states), numerical solution of (14) is 
accomplished with a fourth-order, fixed-step Runge-Kutta 
algorithm. Details of the formulation are given in [6]. 




Figure 3: Inertial Coordinate System and 
Initial Orbit Plane Geometry 


B . Attitude Dynamics 

Input signals for the attitude dynamics model are orbit 
radius vector expressed in inertial coordinates, and aerody- 
namic torque expressed in body coordinates. The body co- 
ordinate system is the same as the geometric coordinate 
system defined in Fig. 2 except that its origin is at the center 
of mass of the spacecraft. Output signals are body rates and 
attitude angles. The mathematical model, of the type found 
in [5], represents rotational motion of a rigid body, and is 
defined by 

/•<u + <5x/-<y = f gg +f aer0 , (16) 

where / is the moment of inertia matrix, and <B is the angu- 
lar velocity vector associated with the body frame. , the 
gravity gradient torque, is given by [7] 

f gg = GTn *f h rXl-r , (17) 

gg ^3 

where r is a unit vector corresponding to R expressed in 
body coordinates. , the aerodynamic torque, is derived 

from output of the aerodynamics subsystem simulation. 
Equation (16) is solved using an Euler angle formulation of 
the kinematics followed by rearrangement to first-order 
form (six states) and numerical integration by Runge-Kutta. 
Attitude angles of Z-Y-X Euler rotation sequence define 
body frame orientation with respect to the inertial frame. 
Details of the formulation are given in [6]. 

C. Aerodynamics 

The aerodynamics simulation input signals are orbit ra- 
dius vector in inertial coordinates and orbit velocity vector 
in body coordinates. Output signals are aerodynamic force 
and torque expressed in geometric coordinates. The space- 
craft is modeled as a box with six rectangular side panels 


(Fig. 2). Drag force for the zth panel is defined empirically 
by [7] 
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where p is atmospheric density, V is velocity magnitude 


relative to the atmosphere, C d is a drag coefficient, A is 

panel area, h is a unit vector normal to the panel and di- 
rected outward from the box, and V is a unit vector corre- 
sponding to V expressed in geometric coordinates. Because 
rotation of the Earth is neglected, no distinction between 
relative and absolute velocity is made. The total aerody- 
namic force is given by 

■ < 19 > 

The aerodynamic torque about the geometric center of 
the box is given by 
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f = (20) 

1=1 

where is a vector directed from the origin of the geo- 
metric coordinate system to the center of pressure of the zth 
panel. Finally, 7^ , as required in (16), is given by 


r aero =r-(r cm xF) , ( 2 i) 

where r cm is a vector directed from the geometric center of 

the spacecraft to its center of mass. Equation (21) is not 
contained within the aerodynamics subsystem simulation, 
but is implemented as a coupler function. 

V. SYSTEM SIMULATION STRUCTURES 

Three system-level simulations of different structure have 
been built. The first simulation, illustrated in Fig. 4, uses an 
implementation of recursive feedback with centralized cou- 
pling, while a second simulation, illustrated in Fig. 5, uses 
distributed coupling. The third simulation uses a conven- 
tional first-order, single-equation-set approach and is the 
standard against which the recursive feedback approaches 
are compared. For the distributively coupled simulation, the 
subsystems have been arranged so that, to the extent possi- 
ble, feedback is minimized. Because there is an information 
loop that ultimately must be closed, it is not possible to 
completely eliminate feedback. The subsimulations are 
identical to those of the centrally coupled system, but the 
system-level numerical process is substantially affected be- 
cause, among other things, the number of scalar signal paths 
that cross the recursive transfer point has changed. 




Recursive Transfer Point *** Spacecraft Geometric Coordinates 

* Inertial Coordinates «■*** From Inertial to Body Frame, 

** Bod y Coordinates (cm origin) 2^ Y -X Rotation Sequence 


Figure 4: Centrally Coupled Simulation Diagram 

Simulations for each of the three subsystems have been 
constructed as functionally separate entities and combined 
within the framework of a single FORTRAN computer pro- 
gram that couples them by way of a recursive feedback 
process. Input and output signals (functions of time) are 
represented by a sequence of linearly connected points 
equally spaced in time, and additional points are obtained 
through interpolation by each simulation, if needed. The 
orbit dynamics and attitude dynamics simulations have 
separate numerical integration processes, while the aerody- 
namics simulation has no integration process because the 
model is purely logical-algebraic. All three simulations are 
“stand-alone” because each, in principle, is capable of pro- 
ducing output signals from input signals and initial condi- 
tions. Initial conditions, of course, do not apply to the aero- 
dynamic simulation because no integration is involved. The 
conventional simulation is implemented from the combined 
set of equations that define the subsystem models. They 
have been collected and rearranged, by hand, to a single set 
of 12 first-order differential equations, an arrangement 
sometimes called a “state variables” approach to system 
formulation. In that context, the state vector is defined as 



Figure 5: Distributively Coupled Simulation Diagram 


X = [fl|V|<f>|£] r , (22) 

where R ,V and £2 are three-element row vectors 
representing orbit radius, velocity, attitude angles, and body 
rates. The system differential equations can now be repre- 
sented as a first-order set by 

X = f(X). (23) 

A fourth-order, fixed-step Runge-Kutta numerical integra- 
tion algorithm specifically designed to solve systems in the 
form of (13) is employed to propagate the system re- 
sponses. 

IV. NUMERICAL RESULTS 

A. Definition of the Run Case 

Spacecraft mass properties, dimensions, and aerody- 
namic data are specified in Fig 2. Initial attitude angles and 
body rates are zero. The orbit is initialized so that in the 
absence of aerodynamic drag it would be circular at an alti- 
tude of 300 km and an inclination of 30 deg. The orbit ini- 
tialization point is in the X-Z plane of the coordinate sys- 
tem shown in Fig. 3. The spacecraft is allowed to tumble 
under the influence of gravity and aerodynamics for a time 
period of 60,000 sec, or about 1 1 orbital revolutions. 

Algorithm parameters in both recursive method simula- 
tions were set to use five sample points per convergence 
interval (including both end points), and the interval length 
was 0.4 sec. Root mean square normalized and absolute 

_1 2 

convergence error tolerances were 10 for the centrally 

coupled case and 10 14 for the distributively coupled case. 
The integration step size for the orbit and attitude dynamics 
sub-simulations was 0.1 sec, and the aerodynamics sub- 
simulation computed aerodynamic effects at a 0.1 -sec sam- 
ple interval. Integration step size for the conventional simu- 
lation was also 0.1 sec. 

B. Comparison of Responses 

Responses from all the simulations agree for about 
45,000 sec of the 60,000-sec duration. After that, the recur- 
sive method responses begin to diverge from the conven- 
tional method. However, the recursive methods continue to 
agree with each other for the full 60,000 sec. Figs. 6-12 
present plots of time segments specifically selected to illus- 
trate the early agreement and subsequent deviation of re- 
sponses. Attitude angles are shown in Figs. 6-8, body 
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Figure 12: Altitude vs. Time 


rates are shown in Figs. 9-11, and altitude is shown for the 
final 1000 sec in Fig. 12. Deviations are attributed primarily 
to the fact that aerodynamic forces and torques are linearly 
interpolated over the subsimulation integration step size in 
the recursive methods, while they are computed at points 
internal to the integration interval in the conventional im- 
plementation. Improvement could likely be found through a 
more sophisticated interpolation method, a shorter conver- 
gence interval, or an increased number of points per con- 
vergence interval for signal definition. It must be remem- 
bered, however, that eventual disagreement of dynamical 
simulations driven by different numerical processes is al- 
ways expected. Additional cases involving two- and three- 
axis attitude control were studied in [6], and no significant 
deviation became apparent for the full 60,000 sec duration. 

C. Performance 

For the example presented, the recursive methods were 
slower than the conventional method by a factor of 4 to 6 
depending mainly on the coupling approach. The centrally 
coupled version was slowest, but no advantage of the pos- 
sibility for parallel execution of subsystem simulations was 
taken in the current single-computer implementation. A 
multiplatform (or multiprocessor) implementation could 
naturally take such advantage; however, that possibility 
does not exist in the distributively coupled version. In the 
recursive methods, response at the latest stage is compared 
with response from the previous stage to determine conver- 
gence with respect to a set of error tolerances. In the con- 
ventional method, no active control of error is present, and 
inclusion of such a mechanism would require additional 
computation. Conversely, predetermination of the number 
of recursion stages and elimination of the convergence 
check could significantly reduce the computation load for 
the recursive methods. While it is recognized that speed is 
important, none of the simulations were refined for effi- 
ciency, and it is believed that many advantages of the re- 
cursive feedback approach lie elsewhere. 


VL SUMMARY AND CONCLUSIONS 

An example nonlinear, continuous simulation of satellite 
motion has been successfully constructed using two varia- 
tions of recursive feedback to couple three separate sub- 
simulations of orbit dynamics, attitude dynamics, and aero- 
dynamics. Results have been verified by direct comparison 
with a conventionally constructed simulation of the same 
system. Each subsimulation deals with one engineering dis- 
cipline, and appropriate interactions are implemented recur- 
sively at the system level. The subsimulations can, in prin- 
ciple, be run separately on remote and dissimilar platforms 
while coupling is implemented by way of network. The sys- 
tem simulation can be designed so that coupling signals 
represent physically identifiable variables associated with 
physical interfaces among the subsystems. The example 
employs three subsimulations; the method framework ac- 
commodates an arbitrary number. Clearly, further investiga- 
tion of recursively coupled subsimulations as a multi- 
discpline system simulation architecture is warranted. 
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